# TODO: create PM image from aod and temp (PM regression)
# 
# Author: phamha
###############################################################################

#import library
library(gstat)
library(RPostgreSQL)
library(stringr)
library(raster)

#import source
setwd("/home/phamha/project/Regression")
source("constant.R")
source("util.R")

data_folder = "/var/www/html/"
met_folder = "/home/phamha/MetTiff/"
myd_folder = "/var/www/html/fimo/apom/Product/MYD01PM/"
mod_folder = "/var/www/html/fimo/apom/Product/MOD01PM/"


mod_query = "SELECT mod04.aqstime, mod04.filename,mod04.filepath,mod07.filename as temp_filename,mod07.filepath as temp_filepath FROM apom.satresampmod04 as mod04 inner join apom.satresampmod07temperature as mod07 ON (mod04.aqstime = mod07.aqstime) where mod04.aqstime < '2014-09-12 00:00:00'::timestamp"
myd_query = "SELECT mod04.aqstime, mod04.filename,mod04.filepath,mod07.filename as temp_filename,mod07.filepath as temp_filepath FROM apom.satresampmyd04 as mod04 inner join apom.satresampmyd07temperature as mod07 ON (mod04.aqstime = mod07.aqstime) where mod04.aqstime < '2014-09-12 00:00:00'::timestamp"

satellite_image_to_pm = function(sate_data,model_type,cook_type,query,out_folder,min_aod,max_aod,min_temp,max_temp){
	data = getDataFromDB(query)
	total_record = nrow(data)
	result_dataframe = data.frame()
	for(i in 1:total_record){
		mod04_aqstime = data$aqstime[i]
		aqstime = strptime(mod04_aqstime,format="%Y-%m-%d %H:%M:%S")
		aqstime = aqstime + 25200
		
		month = format.Date(aqstime,"%m")
		year = format.Date(aqstime,"%Y")
		aod_filename = str_trim(data$filename[i])
		aod_path = str_trim(data$filepath[i])
		aod_path = paste(data_folder,aod_path,sep = "")
		
		point_index = regexpr(".tif",aod_filename) - 1
		file_name = substr(aod_filename,1,point_index)
		aod_mask_file = paste(aod_path,file_name,"_mask.tif",sep = "")
		
		temp_filename = str_trim(data$temp_filename[i])
		temp_path = str_trim(data$temp_filepath[i])
		temp_path = paste(data_folder,temp_path,sep = "")
		
		temp_mask_file = paste(temp_path,temp_filename,"_T_10km_mask.tif",sep = "")
		
		#aod_dataset = raster("/home/phamha/Desktop/aod.tif")
		aod_dataset = raster(aod_mask_file)
		aod_data = values(aod_dataset)
		aod_data = aod_data * 0.00100000004749745
		
		corxy=coordinates(aod_dataset)
		x=corxy[,'x']
		y=corxy[,'y']
		
		#temp_dataset = raster("/home/phamha/Desktop/temp.tif")
		temp_dataset = raster(temp_mask_file)
	
		temp_data = values(temp_dataset)
		temp_data = (temp_data + 15000) * 0.00999999977648258
		
		avg_temp_file  = paste(met_folder,"temp",as.numeric(month),".tif",sep="")
		avg_rh_file    = paste(met_folder,"rh",as.numeric(month),".tif",sep="")
		avg_preci_file = paste(met_folder,"preci",as.numeric(month),".tif",sep="")
		
		
		avg_temp_dataset  =  raster(avg_temp_file)
		avg_rh_dataset    =  raster(avg_rh_file)
		avg_preci_dataset =  raster(avg_preci_file)
		
		avg_temp_data  = values(avg_temp_dataset) 
		avg_rh_data    = values(avg_rh_dataset) 
		avg_preci_data = values(avg_preci_dataset)
		
		
		#select 30%
		#pm_data = regress_predict("mod","linear","4np",aod_data,temp_data,avg_temp_data,avg_rh_data,avg_preci_data)

		pm25 = regress_predict(sate_data,model_type,cook_type,aod_data,temp_data,avg_temp_data,avg_rh_data,avg_preci_data)
		total_pixel = sum(!is.na(aod_data))
		ratio = total_pixel/2024*100
	
		if(ratio>=30){
			table = data.frame(x,y,aod_data,temp_data,avg_temp_data,avg_rh_data,avg_preci_data,pm25)
			valid_table = subset(table,aod_data>=min_aod&aod_data<=max_aod&temp_data>=min_temp&temp_data<=max_temp)
			invalid_table  = subset(table,aod_data<min_aod|aod_data>max_aod|temp_data<min_temp|temp_data>max_temp)
			
			valid_pixel = nrow(valid_table)
			invalid_pixel = nrow(invalid_table)
			result_dataframe = rbind(result_dataframe,data.frame(year,file_name,total_pixel,valid_pixel,invalid_pixel))
			
			table$pm25[table$aod_data<min_aod|table$aod_data>max_aod|table$temp_data<min_temp|table$temp_data>max_temp]<-NA
			og_raster = aod_dataset
			totalCell = ncell(og_raster)
			og_raster[1:totalCell] = table$pm25
			
			out_path = paste(out_folder,year,"/",sep="")
			## valid_table_file = paste(out_path,file_name,"_valid.csv")
			## invalid_table_file = paste(out_path,file_name,"_in_valid.csv")
			pm_file = paste(out_path,file_name,"_og.tif")
			
			## write.csv(valid_table,file=valid_table_file)
			## write.csv(invalid_table,file=invalid_table_file)
			writeRaster(og_raster,filename=pm_file,format="GTiff")
			print(file_name)
			
		}
		
	
	}
	result_file = paste(out_folder,"_range_statis.csv",sep="")
	write.csv(result_dataframe,file=result_file)
}
#satellite_image_to_pm("mod","linear","4np",mod_query,mod_folder,min_aod_mod,max_aod_mod,min_temp_mod,max_temp_mod)
satellite_image_to_pm("myd","linear","4np",myd_query,myd_folder,min_aod_myd,max_aod_myd,min_temp_myd,max_temp_myd)
print("done")
